#include "AnaDef.h"
#include "XGLUtils.h"
#include "AliTRDdEdxBaseUtils.h"
#include "style.h"

TGraphAsymmErrors * GEANTcurve()
{
  /*//from Anton (http://www-alice.gsi.de/trd/papers/dedx/xeco15.r)
    4  xeco15
    1.100000   75.009056
    1.200000   45.508083
    1.300000   35.299252
    1.500000   27.116327
    1.800000   22.734999
    2.000000   21.411915
    2.500000   19.934095
    3.000000   19.449375
    4.000000   19.344431
    7.000000   20.185553
   10.000000   21.027925
   20.000000   22.912676
   40.000000   24.933352
   70.000000   26.504053
  100.000000   27.387468
  300.000000   29.566597
  600.000000   30.353779
 1000.000000   30.787134
 3000.000000   31.129285
10000.000000   31.157350
  */

  TGraphAsymmErrors *gr = new TGraphAsymmErrors;
  gr->SetName("GEANT3");
  gr->SetTitle("PAI model");
  //gr->SetMarkerStyle(34);
  gr->SetMarkerColor(kBlue+1);
  gr->SetLineColor(kBlue+1);
  gr->SetLineStyle(kDotted);
  gr->SetLineWidth(2);

  Int_t ip=0;
  gr->SetPoint(ip++,     1.100000,  75.009056);
  gr->SetPoint(ip++,     1.200000,  45.508083);
  gr->SetPoint(ip++,     1.300000,  35.299252);
  gr->SetPoint(ip++,     1.500000,  27.116327);
  gr->SetPoint(ip++,     1.800000,  22.734999);
  gr->SetPoint(ip++,     2.000000,  21.411915);
  gr->SetPoint(ip++,     2.500000,  19.934095);
  gr->SetPoint(ip++,     3.000000,  19.449375);
  gr->SetPoint(ip++,     4.000000,  19.344431);
  gr->SetPoint(ip++,     7.000000,  20.185553);
  gr->SetPoint(ip++,    10.000000,  21.027925);
  gr->SetPoint(ip++,    20.000000,  22.912676);
  gr->SetPoint(ip++,    40.000000,  24.933352);
  gr->SetPoint(ip++,    70.000000,  26.504053);
  gr->SetPoint(ip++,   100.000000,  27.387468);
  gr->SetPoint(ip++,   300.000000,  29.566597);
  gr->SetPoint(ip++,   600.000000,  30.353779);
  gr->SetPoint(ip++,  1000.000000,  30.787134);
  gr->SetPoint(ip++,  3000.000000,  31.129285);
  gr->SetPoint(ip++, 10000.000000,  31.157350);

  XGLUtils::ScaleGraph(gr, 1/19.344431);

  return gr;
}

TGraphAsymmErrors * Combine1kG5kG(TGraphAsymmErrors * g1, TGraphAsymmErrors * g5)
{
  TGraphAsymmErrors * gout = new TGraphAsymmErrors;

  TGraphAsymmErrors *gtmp = 0x0;
  const Double_t cut=25;
  Int_t ip=0;

  gtmp = g1;
  for(Int_t ii=0; ii<gtmp->GetN(); ii++){
    Double_t xx=-999, yy=-999;
    gtmp->GetPoint(ii,xx,yy);
    
    if(xx>cut)
      break;
    
    printf("Combine1kG5kG %s %d %e %e\n", gtmp->GetName(), ii, xx, yy);

    gout->SetPoint(ip, xx, yy);
    gout->SetPointError(ip, gtmp->GetErrorXlow(ii), gtmp->GetErrorXhigh(ii), gtmp->GetErrorYlow(ii), gtmp->GetErrorYhigh(ii));
    ip++;
  }

  gtmp = g5;
  for(Int_t ii=0; ii<gtmp->GetN(); ii++){
    Double_t xx=-999, yy=-999;
    gtmp->GetPoint(ii,xx,yy);
    
    if(xx<cut)
      continue;
    
    printf("Combine1kG5kG %s %d %e %e\n", gtmp->GetName(), ii, xx, yy);

    gout->SetPoint(ip, xx, yy);
    gout->SetPointError(ip, gtmp->GetErrorXlow(ii), gtmp->GetErrorXhigh(ii), gtmp->GetErrorYlow(ii), gtmp->GetErrorYhigh(ii));
    ip++;
  }

  return gout;
}

void CorrectTestbeam(TGraphAsymmErrors * gr, const Bool_t ktr)
{
  const Int_t nx = gr->GetN();
  Double_t xs[nx], ys[nx], ey[nx], ehy[nx];
  memcpy(xs, gr->GetX(), nx*sizeof(Double_t));
  memcpy(ys, gr->GetY(), nx*sizeof(Double_t));
  memcpy(ey, gr->GetEYlow(), nx*sizeof(Double_t));
  memcpy(ehy, gr->GetEYhigh(), nx*sizeof(Double_t));

  XGLUtils::FFunc func = 0x0;
  Int_t npar = -999;
  Double_t par[10];
  Double_t err[10];
  Double_t chi[10];
  
  if(ktr){
    func = AliTRDdEdxBaseUtils::MeandEdxTR;
    npar = 8;

    Double_t tmp[]={0.706, 1.85, 7.8,0.2, 4.4, 4e-3, 2.26, 0.95};
    memcpy(par, tmp, npar*sizeof(Double_t));
  }
  else{
    func = AliTRDdEdxBaseUtils::MeandEdx;
    npar = 5;

    Double_t tmp[]={0.2, 4.4, 4e-3, 2.26, 0.95};
    memcpy(par, tmp, npar*sizeof(Double_t));
  }

  XGLUtils::ChisquareFit(nx, xs, ys, 0x0, ey, func, npar, par, err, chi);

  for(Int_t ii=0; ii<nx/2; ii++){
    Int_t kk = ii+nx/2;
    if(kk==19){
      ///u/xlu/.task/dedx/testBeam/NewOutputStructure_fitting__Mo_18_Jun_13_10_23_CEST_2012/raw/dEdx_2004_inv8_h009/see.log:AnadEdxFit::Save SPSIONTRh110pdf fitfail 1 0 4
      //the true x[9] is missing, so this x[9] is actually x[10], which should correct x[20]
      kk=20;
    }

    const Double_t scale = func(&(xs[ii]), par)/ys[ii];
    printf("correcting %s %d %e with %d %e -- %e\n", gr->GetName(), kk, xs[kk], ii, xs[ii], scale);
    gr->SetPoint(kk, xs[kk], ys[kk]*scale);
    gr->SetPointError(kk, 0, 0, ey[kk]*scale, ehy[kk]*scale);
  }

  //does not help
  /*
  for(Int_t ii=nx/2; ii<nx; ii++){
    Int_t kk = ii-nx/2;
    if(kk==9){
      ii++;
    }

    const Double_t scale = func(&(xs[ii]), par)/ys[ii];
    printf("correcting %s %d %e with %d %e -- %e\n", gr->GetName(), kk, xs[kk], ii, xs[ii], scale);
    gr->SetPoint(kk, xs[kk], ys[kk]*scale);
    gr->SetPointError(kk, 0, 0, ey[kk]*scale, ehy[kk]*scale);
  }
  */
  printf("\n");
}

void drawSummary()
{
  TList *ll=(TList*)gDirectory->Get("mpv");
  if(!ll){
    printf("no mpv!\n"); gDirectory->ls(); exit(1);
  }

  const TString gname[]={
    "grSPSpureIONmpv",   //0
    "grSPSIONTRmpv",     //1
    "grppmpv",           //2
    "grCosmic5kGINmpv",  //3
    "grCosmic5kGOUTmpv", //4 -> grCosmicOUTmpv
    "grCosmic1kGOUTmpv", //5 -> grCosmic5kGOUTmpv
    "grCosmicOUTmpv"     //6 -> grCosmic1kGOUTmpv
  };      
  const Int_t nraw=sizeof(gname)/sizeof(TString);
  TGraphAsymmErrors *grs[nraw];
  for(Int_t ii=0; ii<=5; ii++){
    grs[ii]=(TGraphAsymmErrors*)ll->FindObject(gname[ii]);
    if(!grs[ii]){
      printf("gr not found %s %d\n", gname[ii].Data(), ii);exit(1);
    }

    if(gname[ii].Contains("SPS")){
      CorrectTestbeam(grs[ii], gname[ii].Contains("TR"));
    }
  }

  TGraphAsymmErrors * tmp4 = Combine1kG5kG(grs[5], grs[4]);
  grs[6] = grs[5];
  grs[5] = grs[4];
  grs[4] = tmp4;

  grs[4]->SetName("grCosmicOUTmpv");
  const Int_t ng = 5;

  //===========================================================
  //===========================================================

  TCanvas *c1=new TCanvas("c1","",700,500);
  PadSetup(c1);
  c1->SetLeftMargin(0.12);
  c1->SetBottomMargin(0.12);
  c1->SetRightMargin(0.05);
  c1->SetTopMargin(0.05);

  gPad->SetLogx();
 
  TLegend *lg = new TLegend(0.55,0.16,0.94,0.52);
  lg->SetName("lg");
  ResetStyle(lg, 0.1);
  //lg->SetBorderSize(1);

  const TString tits[]={
    "#pi, e, dE/dx (testbeam)", 
    "#pi, e, dE/dx+TR (testbeam)", 
    "p, #pi, e (pp #sqrt{s}=7 TeV)",
    "#mu, dE/dx (cosmic rays)",
    "#mu, dE/dx+TR (cosmic rays)"    
  };

  const Int_t sty[] = {20, 21, 29, 26, 32, 27, 24, 30};
  const Int_t col[] = {kBlue, kMagenta, kGreen+3, kRed,  kBlack, kBlack, kYellow, kRed};

  const Double_t xmin = pow(10,-0.5);
  const Double_t xmax = pow(10,4.6);
  Int_t iselected=0;
  for(Int_t ii=0; ii<ng; ii++){
    ResetStyle(grs[ii]);
    grs[ii]->SetTitle(tits[ii]);
    grs[ii]->GetXaxis()->SetTitle("#beta#gamma");
    grs[ii]->GetYaxis()->SetTitle("most probable TRD signal (a.u.)");
    grs[ii]->SetMarkerStyle(sty[ii]);
    grs[ii]->SetMarkerColor(col[ii]);
    grs[ii]->SetLineColor(col[ii]);
    grs[ii]->GetXaxis()->SetLimits(xmin, xmax);
    grs[ii]->GetYaxis()->SetRangeUser(0, 2.8);

    grs[ii]->Draw(iselected?"p":"ap");

    lg->AddEntry(grs[ii], grs[ii]->GetTitle(), "p");

    lg->Draw();

    c1->Print(Form("mpv%d.eps",iselected++));
  }

  TLatex *tt=new TLatex(0.35, 0.88, "ALICE TRD (Xe-CO_{2} [85-15])");
  tt->SetName("tt");
  ResetStyle(tt);
  tt->Draw();
  //=================================================================
  //=================================================================
  
  Int_t ngr = 0;
  Int_t ndata = 0;
  XGLUtils::FFunc func = 0x0;
  Int_t npar = 0;  
  Double_t xs[500], ys[500], ex[500], ey[500];
  Double_t par[10], err[10], chi[10];
  const Double_t fCut1 = 25;
  const Double_t fCut2 = 2.93e3;
  Double_t mins[]={-1e10, fCut1, fCut2};
  Double_t maxs[]={fCut1, fCut2, 1e10};

  //=================================================================
  TGraphAsymmErrors *grdEdx[]={grs[2], grs[3], grs[0]};
  ngr = sizeof(grdEdx)/sizeof(TGraphAsymmErrors *);
  ndata = XGLUtils::GraphToData(ngr, grdEdx, mins, maxs, xs, ys, ex, ey);
  func = AliTRDdEdxBaseUtils::MeandEdx;
  npar = 5;
  Double_t tmp1[]={0.2, 4.4, 4e-3, 2.26, 0.95}; 
  memcpy(par, tmp1, npar*sizeof(Double_t));
  XGLUtils::ChisquareFit(ndata, xs, ys, ex, ey, func, npar, par, err, chi);

  TH1D *h1 = XGLUtils::GetHfit("MeandEdx", func, par, xmin, xmax, kTRUE);
  h1->SetLineColor(kBlack);
  h1->SetLineWidth(1);
  h1->SetLineStyle(kDashed);
  h1->Draw("same hist C");
  lg->AddEntry(h1,"fit, ALEPH parametrization","l");
  c1->Print(Form("mpv%d.eps",iselected++));

  //=================================================================
  TGraphAsymmErrors *grAll[]={grs[2], grs[4], grs[1]};
  ngr = sizeof(grAll)/sizeof(TGraphAsymmErrors *);
  ndata = XGLUtils::GraphToData(ngr, grAll, mins, maxs, xs, ys, ex, ey);
  func = AliTRDdEdxBaseUtils::MeandEdxTR;
  npar = 8;
  Double_t tmp2[]={0.706, 1.85, 7.8, par[0], par[1], par[2], par[3], par[4]};
  //Bool_t pfix[]={0,0,0,1,1,1,1,1};
  memcpy(par, tmp2, npar*sizeof(Double_t));
  XGLUtils::ChisquareFit(ndata, xs, ys, ex, ey, func, npar, par, err, chi);//, pfix);

  TH1D *h2 = XGLUtils::GetHfit("MeandEdxTR", func, par, xmin, xmax, kTRUE);
  h2->SetLineColor(kRed);
  h2->SetLineWidth(1);
  h2->SetLineStyle(kSolid);
  h2->Draw("same hist C");
  lg->AddEntry(h2,"fit, ALEPH param. + logistic f.","l");
  c1->Print(Form("mpv%d.eps",iselected++));

  //=================================================================
  DrawLogo(c1, 0.135, 0.21, "19/08/2011");
  c1->Print(Form("mpv_logo.eps"));
  
  //=================================================================
  /*
  TGraphAsymmErrors * geant = GEANTcurve();
  geant->Draw("l");
  lg->AddEntry(geant,"PAI model (Geant3)","l");
  c1->Print(Form("mpv%d.eps",iselected++));
  */
  
  TFile *fout=new TFile("mpv.root", "recreate");
  for(Int_t ig=0; ig<ng; ig++){
    grs[ig]->Write();
  }
  //geant->Write();
  tt->Write();
  lg->Write();
  h1->Write();
  h2->Write();
  c1->Write();
  fout->Save();
  fout->Close();

  //delete geant;
  delete tt;
  delete lg;
  delete h1;
  delete h2;
  delete c1;
}

int main(int argc, char * argv[])
{
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  SetGlobalStyle(0,1);

  TFile::Open(argv[1]);
  drawSummary();

  return 0;
}
